library(tidyverse)
library(sf)
library(exactextractr)
library(patchwork)
library(ggtext)
library(units)
library(modelsummary)
library(gt)
library(rayshader)
set.seed(123)
distrito <- read_sf("dados/distrito/SIRGAS_SHP_distrito.shp") %>% st_set_crs("epsg:31983")
Microdados do Censo de 2022
# Read dados do censo 2022
censo <- read_sf("dados/censo/SP_Malha_Preliminar_2022.shp") %>%
filter(CD_MUN == "3550308") %>%
st_transform("epsg:31983") %>% # Sistema de coordenadas do geosampa
select(id_setor_censitario = CD_SETOR, v0001:v0007) %>%
mutate(area_setor = st_area(geometry) %>% as.numeric())
censo %>%
st_drop_geometry() %>%
modelsummary::datasummary_skim()
tinytable_prpgm970f1dbp778w210
| |
Unique |
Missing Pct. |
Mean |
SD |
Min |
Median |
Max |
Histogram |
| v0001 |
1225 |
0 |
415.0 |
230.9 |
0.0 |
397.0 |
2221.0 |
 |
| v0002 |
592 |
0 |
181.1 |
94.3 |
0.0 |
177.0 |
1372.0 |
 |
| v0003 |
588 |
0 |
180.9 |
94.3 |
0.0 |
177.0 |
1371.0 |
 |
| v0004 |
26 |
0 |
0.2 |
1.6 |
0.0 |
0.0 |
170.0 |
 |
| v0005 |
14492 |
0 |
2.6 |
0.6 |
0.0 |
2.7 |
10.3 |
 |
| v0006 |
7559 |
0 |
13.6 |
12.1 |
0.0 |
10.8 |
100.0 |
 |
| v0007 |
511 |
0 |
156.4 |
81.9 |
0.0 |
152.0 |
950.0 |
 |
| area_setor |
27592 |
0 |
55126.5 |
360317.6 |
196.3 |
23142.2 |
22848344.5 |
 |
gg <- ggplot() +
geom_sf(data = censo %>%
mutate(densidade = v0001/area_setor,
decil = ntile(densidade, 10)),
aes(fill = factor(decil)), color = NA) +
scale_fill_viridis_d() +
labs(fill = "Decil de densidade") +
theme_void() +
theme(legend.position = c(.8, .3))
ggsave("tex/imagens/mapa.pdf", gg, width = 7, height = 12)
gg

censo %>%
st_drop_geometry() %>%
summarize("Total de Pessoas" = sum(v0001),
"Total de Domicílios" = sum(v0002),
"Total de Domicílios Particulares Ocupados" = sum(v0007)) %>%
pivot_longer(everything(), names_to = "Variável", values_to = "Valor") %>%
gt()
| Variável |
Valor |
| Total de Pessoas |
11451999 |
| Total de Domicílios |
4996529 |
| Total de Domicílios Particulares Ocupados |
4316336 |
gg <- censo %>%
mutate(distancia_centro = st_distance(geometry,
read_sf("dados/distrito/SIRGAS_SHP_distrito.shp") %>%
st_set_crs("epsg:31983") %>%
filter(ds_nome == "SE")) %>%
as.numeric() %>%
cut(breaks = 10^3*(0:40), labels = FALSE)) %>%
st_drop_geometry() %>%
group_by(distancia_centro) %>%
summarize(densidade = sum(v0001) / sum(as.numeric(area_setor) / 10^6)) %>%
mutate(teorico = 2*10^4*exp(-distancia_centro/10)) %>%
ggplot(aes(x = distancia_centro)) +
geom_col(aes(y = densidade), fill = "#2BAA92FF") +
geom_line(aes(y = teorico), lwd = .8, linetype = "dashed", color = "#D32934FF") +
scale_y_continuous(labels = scales::comma_format(big.mark = ".")) +
scale_x_continuous(labels = scales::comma_format(big.mark = ".")) +
labs(x = "Distância ao centro em km",
y = expression("Densidade populacional por" ~ km^2)) +
theme_classic()
ggsave("tex/imagens/densidade_distcentro.pdf", gg, width = 5, height = 3)
gg

Base de Dados do IPTU
IPTU <- read.csv("dados/IPTU/IPTU_2024.csv", sep=";", encoding = "latin1") %>%
as_tibble() %>%
select(sql = "NUMERO.DO.CONTRIBUINTE",
condominio = "NUMERO.DO.CONDOMINIO",
area_terreno = "AREA.DO.TERRENO",
area_construida = "AREA.CONSTRUIDA",
area_ocupada = "AREA.OCUPADA",
pavimentos = "QUANTIDADE.DE.PAVIMENTOS",
ano_construcao = "ANO.DA.CONSTRUCAO.CORRIGIDO",
tipo = "TIPO.DE.PADRAO.DA.CONSTRUCAO") %>%
# Separação do número de contribuinte (SQL) em setor quadra e lote
mutate(setor = str_sub(sql, 1, 3),
quadra = str_sub(sql, 4, 6),
# Quando o lote é um condomínio, haverá vários SQLs no mesmo lote. CD = Condomínio
lote = str_sub(sql, 7, 10) %>%
ifelse(condominio == "00-0", ., paste("CD", str_sub(condominio, 1, 2), sep = "")),
# Tipo de uso
residencial = str_detect(tipo, "Residencial"))
IPTU %>%
select(-c(condominio, tipo)) |>
modelsummary::datasummary_skim()
tinytable_99qstmk1miicw6pgvwz3
| |
Unique |
Missing Pct. |
Mean |
SD |
Min |
Median |
Max |
Histogram |
| area_terreno |
9304 |
0 |
3546.7 |
15146.5 |
1.0 |
800.0 |
4307493.0 |
 |
| area_construida |
6718 |
0 |
157.1 |
1068.9 |
0.0 |
102.0 |
950328.0 |
 |
| area_ocupada |
5159 |
0 |
1276.6 |
3170.0 |
0.0 |
378.0 |
320000.0 |
 |
| pavimentos |
52 |
0 |
9.4 |
8.7 |
0.0 |
5.0 |
98.0 |
 |
| ano_construcao |
127 |
0 |
1948.3 |
283.9 |
0.0 |
1987.0 |
2023.0 |
 |
| residencial |
N |
% |
|
|
|
|
|
|
| FALSE |
455084 |
14.7 |
|
|
|
|
|
|
| TRUE |
2641635 |
85.3 |
|
|
|
|
|
|
gg.right <- IPTU %>%
mutate(uso = case_when(str_detect(tipo, "Residencial") ~ "Residencial",
str_detect(tipo, "Comercial") ~ "Comercial",
str_detect(tipo, "Oficina") ~ "Industria",
str_detect(tipo, "TERRENO") ~ "Terreno",
str_detect(tipo, "Clube") ~ "Entretenimento",
.default = "Outros") %>% as.factor(),
padrao = case_when(str_detect(tipo, "A$") ~ "A",
str_detect(tipo, "B$") ~ "B",
str_detect(tipo, "C$") ~ "C",
str_detect(tipo, "D$") ~ "D",
str_detect(tipo, "E$") ~ "E",
.default = "NA"),
nome = paste("(", padrao, ") ", uso, sep = "")) %>%
group_by(uso, padrao, nome) %>%
summarize(area = sum(area_construida)) %>%
ungroup() %>%
mutate(percentual = area * 100/ sum(area, na.rm = TRUE),
texto = case_when(percentual < 1 ~ "<1%",
percentual < 5 ~ paste(round(percentual,1), "%", sep = ""),
.default = paste(round(percentual, 2), "%", sep = "")),
color = ifelse(uso == "Outros", "white", "black")) %>%
ggplot() +
treemapify::geom_treemap(aes(fill = uso, area = area), size = 2, color = NA) +
treemapify::geom_treemap(aes(alpha = padrao, area = area), fill = "black", color = NA) +
treemapify::geom_treemap_text(aes(area = area, label = nome, color = color), alpha = .4, grow = FALSE, size = 8) +
treemapify::geom_treemap_text(aes(area = area, label = texto, color = color),
alpha = .4, place = "middle", size = 10) +
scale_fill_manual(values = c("Residencial" = "#FAE48B",
"Comercial" = "#A6EEE6",
"Industria" = "#84BD00",
"Entretenimento" = "#FB6467",
"Outros" = "#1A5354")) +
scale_alpha_manual(values = c("A" = 0, "B" = .05, "C" = .1, "D" = .15, "E" = .2, "NA" = 0)) +
scale_color_manual(values = c("black" = "black", "white" = "white")) +
labs(fill = "Tipo de uso", alpha = "Padrão de uso") +
theme(aspect.ratio=1)
gg.left <- IPTU %>%
mutate(uso = case_when(str_detect(tipo, "Residencial") ~ "Residencial",
str_detect(tipo, "Comercial") ~ "Comercial",
str_detect(tipo, "Oficina") ~ "Industria",
str_detect(tipo, "TERRENO") ~ "Terreno",
str_detect(tipo, "Clube") ~ "Entretenimento",
.default = "Outros") %>% as.factor()) %>%
group_by(uso) %>%
summarize(area = sum(area_construida, na.rm = TRUE)) %>%
ungroup() %>%
arrange(desc(uso)) %>%
mutate(percentual = area * 100 / sum(area, na.rm = TRUE),
texto = ifelse(percentual > 5, paste(round(percentual, 1), "%", sep = ""), ""),
ytext = (cumsum(area) + lag(cumsum(area)))/ 2) %>%
ggplot() +
geom_col(aes(y = area, fill = uso, x = "")) +
geom_text(aes(y = ytext, fill = uso, x = "", label = texto), alpha = .5) +
scale_fill_manual(values = c("Residencial" = "#FAE48B",
"Comercial" = "#A6EEE6",
"Industria" = "#84BD00",
"Entretenimento" = "#FB6467",
"Outros" = "#1A5354")) +
theme_void() +
theme(legend.position = "none")
ggsave("tex/imagens/tree_area_construida.pdf",
(gg.left | gg.right) +
plot_layout(widths = c(1,6)),
width = 6.75, height = 5)
(gg.left | gg.right)

IPTU <- IPTU %>%
# Agregar por SQL
group_by(setor, quadra, lote) %>%
summarize(unidades = n(),
area_terreno = median(area_terreno),
area_construida = sum(area_construida),
area_ocupada = median(area_ocupada),
pavimentos = median(pavimentos),
ano_construcao = median(ano_construcao),
residencial = median(residencial)) %>%
ungroup() %>%
mutate(CA = area_construida / area_terreno,
cota_parte = area_terreno / unidades,
verticalizacao = area_construida / area_ocupada) %>%
filter(area_ocupada > 0)
IPTU %>%
filter(residencial == 1) %>%
modelsummary::datasummary_skim()
tinytable_vyqdvjl3txhlyt278vut
| |
Unique |
Missing Pct. |
Mean |
SD |
Min |
Median |
Max |
Histogram |
| unidades |
545 |
0 |
2.4 |
16.4 |
1.0 |
1.0 |
2861.0 |
 |
| area_terreno |
5736 |
0 |
248.8 |
1917.0 |
6.0 |
160.0 |
1442400.0 |
 |
| area_construida |
12097 |
0 |
306.3 |
1954.9 |
2.0 |
120.0 |
400000.0 |
 |
| area_ocupada |
3403 |
0 |
118.4 |
407.8 |
1.0 |
90.0 |
200000.0 |
 |
| pavimentos |
51 |
0 |
1.8 |
1.8 |
1.0 |
2.0 |
52.0 |
 |
| ano_construcao |
166 |
0 |
1981.8 |
14.6 |
1800.0 |
1981.0 |
2023.0 |
 |
| residencial |
1 |
0 |
1.0 |
0.0 |
1.0 |
1.0 |
1.0 |
 |
| CA |
100483 |
0 |
0.9 |
0.8 |
0.0 |
0.8 |
36.1 |
 |
| cota_parte |
17359 |
0 |
208.3 |
1624.9 |
1.5 |
155.0 |
1442400.0 |
 |
| verticalizacao |
53664 |
0 |
1.7 |
1.7 |
0.0 |
1.4 |
180.0 |
 |
IPTU %>%
mutate(Tipo = ifelse(unidades == 1, "Unidade", "Condomínio")) %>%
group_by(Tipo) %>%
summarize("Lotes" = n(),
"Unidades" = sum(unidades)) %>%
gt()
| Tipo |
Lotes |
Unidades |
| Condomínio |
33116 |
1781971 |
| Unidade |
1255235 |
1255235 |
gg <- IPTU %>%
filter(residencial == 1) %>%
filter(cota_parte < 600, CA < 5, verticalizacao < 6) %>%
select("Densidade habitacional" = cota_parte,
"Densidade Construtiva" = CA,
"Pavimentos" = verticalizacao) %>%
pivot_longer(everything()) %>%
ggplot() +
geom_histogram(aes(x = value)) +
facet_wrap(~ name, scales = "free") +
labs(x = "", y = "Frequência") +
theme_classic()
ggsave("tex/imagens/indicadores.pdf", gg, width = 9, height = 3)
gg

gg <- IPTU %>%
filter(residencial == 1) %>%
ggplot(aes(x = verticalizacao, y = CA)) +
geom_bin2d(binwidth = c(1,1) * .5) +
geom_abline(intercept = 0, slope = 1, color = "black", lwd = 0.8, linetype = "dashed") +
geom_smooth(aes(color = "Associação")) +
theme_classic() +
scale_color_manual(values = c("Associação" = "red")) +
scale_fill_gradient(high = "darkblue", low = "white", trans = "log10",
labels = scales::label_number(), breaks = 10^c(2, 3, 4, 5)) +
coord_cartesian(xlim = c(0,15), ylim = c(0,15)) +
labs(fill = "Número de ocorrências", colour = "", x = "Verticalização",
y = "Densidade Construtiva (CA observado)") +
theme(aspect.ratio = 1)
ggsave("tex/imagens/ca_vs_verticalizacao.pdf", gg, width = 7.5, height = 5)
gg

gg <- IPTU %>%
filter(residencial == 1, !is.na(ano_construcao), ano_construcao > 1950) %>%
mutate(ano_construcao = ano_construcao %>% round()) %>%
group_by(ano_construcao) %>%
# Cria quantis, o primeiro e último são o mínimo e máximo, respectivamente
summarize(cota_parte = quantile(cota_parte, 0:6/6),
CA = quantile(CA, 0:6/6),
verticalizacao = quantile(verticalizacao, 0:6/6)) %>%
mutate(quantil = row_number() - 1,
ano_construcao = as.Date(as.character(ano_construcao), format = "%Y")) %>%
pivot_longer(cota_parte:verticalizacao) %>%
group_by(ano_construcao, name) %>%
mutate(ymin = lag(value) %>% replace_na(0)) %>%
filter(quantil %in% 2:5) %>%
ggplot() +
geom_ribbon(aes(x = ano_construcao, ymax = value, ymin = ymin, fill = factor(quantil)), color = "black") +
facet_grid(rows = "name", scales = "free_y",
labeller = as_labeller(c("CA" = "Densidade Construtiva",
"verticalizacao" = "Pavimentos",
"cota_parte" = "Densidade Habitacional"))) +
scale_x_date(limits = as.Date(c("1960", "2022"), format = "%Y")) +
scale_fill_viridis_d(labels = c("20%", "40%", "60%", "80%")) +
theme_bw() +
labs(fill = "Quantil", x = "Ano da construção")
ggsave("tex/imagens/indicadores_tempo.pdf", gg, width = 8, height = 9)
ggsave("tex/imagens/indicadores_tempo_small.pdf", gg, width = 5.5, height = 6)
gg

Geosampa
# Extração do zip file com lotes de cada bairro
if (!"zip" %in% list.files(path = "dados/lotes")){
for (file in list.files(path="dados/lotes/zip", full.names = FALSE) %>%
str_remove("\\.zip")){
unzip(paste("dados/lotes/zip/", file, ".zip", sep = ""),
paste(file, "/", file, ".shp", sep = ""),
exdir = "dados/lotes/unzip")
unzip(paste("dados/lotes/zip/", file, ".zip", sep = ""),
paste(file, "/", file, ".dbf", sep = ""),
exdir = "dados/lotes/unzip")
unzip(paste("dados/lotes/zip/", file, ".zip", sep = ""),
paste(file, "/", file, ".shx", sep = ""),
exdir = "dados/lotes/unzip")
}
}
# Junção dos shapes dos lotes de cada bairro em uma tabela
lotes <- list.files(path="dados/lotes/unzip", full.names = FALSE) %>%
paste("dados/lotes/unzip/", ., "/", ., ".shp", sep = "") %>%
lapply(read_sf) %>%
bind_rows %>%
st_set_crs("epsg:31983")
lotes %>%
mutate(condominio = case_when(lo_condomi == "00" ~ "Unidade",
lo_condomi != "00" ~ "Condomínio"),
tipo = case_when(lo_tp_lote == "F" ~ "Fiscal",
lo_tp_lote == "M" ~ "Espaço livre",
lo_tp_lote == "V" ~ "Via de acesso",
.default = NA),
area = st_area(geometry) %>% as.numeric()) %>%
st_drop_geometry() %>%
ggplot() +
geom_violin(aes(x = factor(condominio), y = area, fill = tipo)) +
scale_y_log10(labels = scales::comma_format(big.mark = ".")) +
labs(title = "Área dos lotes em SP" , x = "", y = "Área em metros quadrados", fill = "Tipo de lote") +
theme_classic()

lotes <- lotes %>%
filter(lo_tp_lote == "F") %>% # Seleção apenas de lotes fiscais
mutate(lo_lote = ifelse(lo_lote == "0000", paste("CD", lo_condomi, sep = ""), lo_lote)) %>%
select(setor = lo_setor, quadra = lo_quadra, lote = lo_lote)
quadras <- read_sf("dados/quadras/SIRGAS_SHP_quadraMDSF.shp") %>%
st_set_crs("epsg:31983") %>%
filter(qd_tipo == "F") %>% # Apenas quadras fiscais
select(setor = qd_setor, quadra = qd_fiscal) %>%
group_by(setor, quadra) %>%
summarize(geometry = st_union(geometry)) %>%
ungroup()
setores <- read_sf("dados/setor/SIRGAS_SHP_setorfiscal.shp") %>%
st_set_crs("epsg:31983") %>%
select(setor = st_codigo) %>%
group_by(setor) %>%
summarize(geometry = st_union(geometry)) %>%
ungroup()
Join IPTU - GeoSampa
# Join dos lotes com IPTU com base no SQL
IPTU.lote <- IPTU %>%
filter(residencial == 1) %>%
left_join(lotes, by = join_by(setor, quadra, lote)) %>%
ungroup()
IPTU.lote %>%
modelsummary::datasummary_skim()
tinytable_c62mbnjfw0fyeogltc0b
| |
Unique |
Missing Pct. |
Mean |
SD |
Min |
Median |
Max |
Histogram |
| unidades |
545 |
0 |
2.4 |
16.4 |
1.0 |
1.0 |
2861.0 |
 |
| area_terreno |
5736 |
0 |
248.8 |
1917.0 |
6.0 |
160.0 |
1442400.0 |
 |
| area_construida |
12097 |
0 |
306.3 |
1954.9 |
2.0 |
120.0 |
400000.0 |
 |
| area_ocupada |
3403 |
0 |
118.4 |
407.8 |
1.0 |
90.0 |
200000.0 |
 |
| pavimentos |
51 |
0 |
1.8 |
1.8 |
1.0 |
2.0 |
52.0 |
 |
| ano_construcao |
166 |
0 |
1981.8 |
14.6 |
1800.0 |
1981.0 |
2023.0 |
 |
| residencial |
1 |
0 |
1.0 |
0.0 |
1.0 |
1.0 |
1.0 |
 |
| CA |
100483 |
0 |
0.9 |
0.8 |
0.0 |
0.8 |
36.1 |
 |
| cota_parte |
17359 |
0 |
208.3 |
1624.9 |
1.5 |
155.0 |
1442400.0 |
 |
| verticalizacao |
53664 |
0 |
1.7 |
1.7 |
0.0 |
1.4 |
180.0 |
 |
IPTU.quadra <- IPTU %>%
filter(residencial == 1) %>%
group_by(setor, quadra) %>%
summarize(unidades = sum(unidades),
area_construida = sum(area_construida),
verticalizacao = weighted.mean(verticalizacao, area_terreno),
area_terreno = sum(area_terreno)) %>%
left_join(quadras, by = join_by(setor, quadra)) %>%
ungroup()
IPTU.quadra %>%
modelsummary::datasummary_skim()
tinytable_9o9p6guvmy33zmw6h0y6
| |
Unique |
Missing Pct. |
Mean |
SD |
Min |
Median |
Max |
Histogram |
| unidades |
976 |
0 |
77.6 |
160.6 |
1.0 |
34.0 |
4762.0 |
 |
| area_construida |
15166 |
0 |
9976.7 |
19124.6 |
8.0 |
4813.0 |
521949.0 |
 |
| verticalizacao |
33567 |
0 |
2.4 |
3.0 |
0.9 |
1.5 |
81.6 |
 |
| area_terreno |
14846 |
0 |
8104.0 |
13914.7 |
50.0 |
6029.0 |
1447291.0 |
 |
IPTU.setor <- IPTU %>%
filter(residencial == 1) %>%
group_by(setor) %>%
summarize(unidades = sum(unidades),
area_construida = sum(area_construida),
verticalizacao = weighted.mean(verticalizacao, area_terreno),
area_terreno = sum(area_terreno)) %>%
left_join(setores, by = join_by(setor)) %>%
ungroup()
IPTU.setor %>%
modelsummary::datasummary_skim()
tinytable_twn6y3cdui0u5nyy0h8d
| |
Unique |
Missing Pct. |
Mean |
SD |
Min |
Median |
Max |
Histogram |
| unidades |
168 |
0 |
15827.8 |
9045.2 |
133.0 |
13694.0 |
46511.0 |
 |
| area_construida |
168 |
0 |
2034480.2 |
1287662.5 |
16480.0 |
1682702.5 |
6822366.0 |
 |
| verticalizacao |
168 |
0 |
3.5 |
2.2 |
1.1 |
2.7 |
12.5 |
 |
| area_terreno |
168 |
0 |
1652592.8 |
1064452.3 |
30130.0 |
1410261.5 |
5087896.0 |
 |
list(IPTU %>% filter(residencial == 1) %>% mutate(base = "bruta"),
IPTU.setor %>% mutate(base = "Setor") %>% filter(!geometry %>% st_is_empty()),
IPTU.quadra %>% mutate(base = "Quadra") %>% filter(!geometry %>% st_is_empty()),
IPTU.lote %>% mutate(base = "Lote") %>% filter(!geometry %>% st_is_empty())) %>%
lapply(function(x) x %>%
group_by(base) %>%
summarize(unidades = sum(unidades))) %>%
bind_rows %>%
pivot_wider(names_from = base, values_from = unidades) %>%
pivot_longer(2:4) %>%
mutate(missing = bruta - value,
missing_percent = (100 *missing / bruta) %>% round(2) %>% paste("%")) %>%
select("Critério de join" = name,
"Erros (unidades)" = missing,
"Percentual de erro" = missing_percent) %>%
gt()
| Critério de join |
Erros (unidades) |
Percentual de erro |
| Setor |
0 |
0 % |
| Quadra |
820 |
0.03 % |
| Lote |
44319 |
1.67 % |
gg <- IPTU.lote %>%
mutate(distancia_centro = st_distance(geometry,
read_sf("dados/distrito/SIRGAS_SHP_distrito.shp") %>%
st_set_crs("epsg:31983") %>%
filter(ds_nome == "SE")) %>%
as.numeric() %>%
cut(breaks = 10^3*(0:40), labels = FALSE)) %>%
st_drop_geometry() %>%
group_by(distancia_centro) %>%
summarize(verticalizacao = weighted.mean(verticalizacao, area_terreno)) %>%
ggplot(aes(x = distancia_centro, y = verticalizacao)) +
geom_col() +
geom_hline(yintercept = 1, linetype = "dotted") +
labs(x = "Distância ao centro em quilômetros",
y = "Verticalização") +
scale_y_continuous(breaks = (1:6)) +
theme_classic()
ggsave("tex/imagens/verticalizacao.pdf", gg, width = 6, height = 4)
gg

gg.lotes <- IPTU.quadra %>%
ggplot() +
geom_sf(data = distrito, fill = "black") +
geom_sf(aes(fill = verticalizacao, geometry = geometry), color = NA) +
theme_void() +
scale_fill_gradient(low = "#004D40FF", high = "#E0F2F1FF")
gg.setor <- IPTU.setor %>%
ggplot() +
geom_sf(data = distrito, fill = "black") +
geom_sf(aes(fill = verticalizacao, geometry = geometry), color = NA) +
theme_void() +
scale_fill_gradient(low = "#004D40FF", high = "#E0F2F1FF")
ggsave("tex/imagens/mapa_verticalizacao.pdf",
(gg.lotes | gg.setor), width = 20, height = 20)
(gg.lotes | gg.setor)

Join IPTU - Censo
# Join dados do IPTU com do Censo através da intersecção das geometrias
IPTU.censo <- censo %>%
st_intersection(IPTU.lote %>% st_as_sf()) %>%
as_tibble() %>%
rename(geometria_intersec = geometry) %>%
# Retomada das geometrias do setor e lotes
left_join(censo %>%
as_tibble() %>%
select(id_setor_censitario,
geometria_setor_censitario = geometry),
by = join_by(id_setor_censitario)) %>%
left_join(IPTU.lote %>%
as_tibble() %>%
select(setor, quadra, lote,
geometria_lote = geometry),
by = join_by(setor, quadra, lote)) %>%
# Cálculo de quanto % do lote está dentro do setor
mutate(percent_intersec = as.numeric(st_area(geometria_intersec) / st_area(geometria_lote)))
IPTU.censo %>%
modelsummary::datasummary_skim()
tinytable_qu3wg7ju8hn9ahiiaw5g
| |
Unique |
Missing Pct. |
Mean |
SD |
Min |
Median |
Max |
Histogram |
| v0001 |
1120 |
0 |
524.4 |
202.6 |
0.0 |
511.0 |
2221.0 |
 |
| v0002 |
543 |
0 |
227.8 |
87.1 |
0.0 |
222.0 |
1372.0 |
 |
| v0003 |
538 |
0 |
227.5 |
87.1 |
0.0 |
221.0 |
1371.0 |
 |
| v0004 |
24 |
0 |
0.2 |
1.0 |
0.0 |
0.0 |
170.0 |
 |
| v0005 |
11463 |
0 |
2.7 |
0.3 |
0.0 |
2.7 |
8.0 |
 |
| v0006 |
6421 |
0 |
13.4 |
10.3 |
0.0 |
11.2 |
100.0 |
 |
| v0007 |
471 |
0 |
194.7 |
72.7 |
0.0 |
190.0 |
950.0 |
 |
| area_setor |
18531 |
0 |
52804.2 |
74529.5 |
196.3 |
41868.0 |
6945494.4 |
 |
| unidades |
541 |
0 |
3.7 |
33.0 |
1.0 |
1.0 |
2861.0 |
 |
| area_terreno |
5624 |
0 |
326.2 |
3606.1 |
6.0 |
162.0 |
1442400.0 |
 |
| area_construida |
11985 |
0 |
460.5 |
4038.2 |
2.0 |
120.0 |
400000.0 |
 |
| area_ocupada |
3350 |
0 |
140.3 |
1052.4 |
1.0 |
90.0 |
200000.0 |
 |
| pavimentos |
51 |
0 |
1.9 |
2.2 |
1.0 |
2.0 |
52.0 |
 |
| ano_construcao |
165 |
0 |
1981.7 |
14.6 |
1800.0 |
1981.0 |
2023.0 |
 |
| residencial |
1 |
0 |
1.0 |
0.0 |
1.0 |
1.0 |
1.0 |
 |
| CA |
99679 |
0 |
0.9 |
0.9 |
0.0 |
0.8 |
36.1 |
 |
| cota_parte |
17099 |
0 |
240.4 |
2779.6 |
1.5 |
156.0 |
1442400.0 |
 |
| verticalizacao |
53205 |
0 |
1.7 |
2.0 |
0.1 |
1.4 |
180.0 |
 |
| percent_intersec |
285710 |
0 |
0.9 |
0.2 |
0.0 |
1.0 |
96.6 |
 |
geometria <- IPTU.censo %>%
filter(id_setor_censitario == "355030810000143P") %>%
select(geometria_setor_censitario) %>%
st_as_sf()
zoom <- 60
bbox <- geometria %>%
st_transform("epsg:31983") %>%
st_bbox()
x_min <- (bbox$xmin - zoom) %>% as.numeric()
y_min <- (bbox$ymin - zoom) %>% as.numeric()
x_max <- (bbox$xmax + zoom) %>% as.numeric()
y_max <- (bbox$ymax + zoom) %>% as.numeric()
corte <- list(xmin = x_min - zoom, ymin = y_min - zoom,
xmax = x_max + zoom, ymax = y_max + zoom)
df <- IPTU.censo %>%
st_as_sf() %>%
st_transform("epsg:31983") %>%
st_crop(corte %>% unlist()) %>%
mutate(percent_intersec = round(percent_intersec * 100))
gg <- df %>%
mutate(texto = ifelse(percent_intersec != 100, percent_intersec, NA),
color = ifelse(percent_intersec == 100, "blue", "red")) %>%
ggplot() +
geom_sf(aes(geometry = geometria_setor_censitario), color = "black", fill = "grey",
lwd = 1.5, linetype ="solid") +
geom_sf(aes(geometry = geometria_intersec, fill = color, color = color),
alpha = .25, linetype = "dashed", lwd = 1) +
geom_sf_label(aes(geometry = geometria_intersec, label = texto), label.size = .01) +
scale_fill_manual(values = c("red" = "red", "blue" = "blue")) +
scale_color_manual(values = c("red" = "red", "blue" = "blue")) +
labs(fill = "Decil de densidade") +
theme_void() +
theme(legend.position = "none") +
coord_sf(xlim = c(x_min, x_max), ylim = c(y_min, y_max))
ggsave("tex/imagens/corte_lote.pdf", gg, width = abs(x_min - x_max)/40, height = abs(y_min- y_max)/40)
gg

Teoricamente o número de domicílios deveria ser próximo do número de
unidades habitacionais enunciadas pelo IPTU. Caso haja mais unidades
habitacionais do que domicílios, houve algum equívoco na medição do
censo, visto que entre domicílios não contam apenas os ocupados. Caso
haja mais domicílios, há unidades que não são contribuintes do IPTU. O
que é preocupante é que casos em que o número é diferente não são
exceções.
gg <- censo %>%
st_drop_geometry() %>%
select(id_setor_censitario, domicilios = v0002) %>%
left_join(IPTU.censo %>%
st_drop_geometry() %>%
filter(residencial == 1) %>%
group_by(id_setor_censitario) %>%
summarize(unidades = sum(unidades * percent_intersec))) %>%
mutate(unidades = replace_na(unidades, 0)) %>%
mutate(espectro_irregularidade = unidades / (unidades + domicilios)) %>%
ggplot() +
geom_histogram(aes(x = espectro_irregularidade)) +
labs(x = "Espectro de irregularidade: unidades / (unidades + domicilios)", y = "Frequência") +
theme_classic() +
theme(plot.caption = element_text(hjust = 0)) +
scale_x_continuous(labels = scales::percent)
ggsave("tex/imagens/disparidade_censoIPTU.pdf", gg, width = 7, height = 3)
gg

Alguns setores censitários não encontram pares na base de lotes. Isso
acontece principalmente por conta de loteamentos irregulares e favelas,
que não são contribuintes do IPTU, e, portanto, não constam na base.
Isso não prejudica a análise, pois estes loteamentos não dependem da
regulamentação urbana, então mudanças nos instrumentos e indicadores não
impactariam essas regiões.
Outras falhas decorrem do erro do join da base do IPTU com lotes,
como apontado anteriormente, mas estes casos são negligenciáveis.
distrito <- read_sf("dados/distrito/SIRGAS_SHP_distrito.shp") %>% st_set_crs("epsg:31983")
lotes_irregulares <- read_sf("dados/lotes_irregulares/SIRGAS_SHP_loteamento.shp") %>%
st_set_crs("epsg:31983")
favela <- read_sf("dados/favela/SIRGAS_SHP_favela.shp") %>%
st_set_crs("epsg:31983")
censo.pontos <- censo %>%
select(id_setor_censitario, geometry, pessoas = v0001) %>%
# Um ponto a cada 1.000 pessoas
mutate(pontos = pessoas %/% 1000 + (runif(length(pessoas)) * 1000 < pessoas %% 1000)) %>%
select(-pessoas) %>%
as_tibble() %>%
left_join(
# Seleção dos setores censitários que não apresentam nenhum contribuinte do IPTU
censo %>%
anti_join(IPTU.censo) %>%
st_drop_geometry() %>%
select(id_setor_censitario) %>%
mutate(erro = TRUE)) %>%
mutate(erro = replace_na(erro, FALSE)) %>%
filter(pontos > 0) %>%
mutate(samples = purrr::map2(geometry, pontos, ~st_sample(.x, size = .y))) %>%
unnest(cols = samples) %>%
as_tibble()
gg <- censo.pontos %>%
ggplot() +
geom_sf(data = distrito, color = NA) +
geom_sf(data = st_union(favela %>% st_union() %>% st_buffer(10),
lotes_irregulares %>% st_union() %>% st_buffer(10)),
aes(fill = "Favelas e lotes irregulares"),
color = "black", size = .1, alpha = .5) +
geom_sf(aes(geometry = samples, color = erro), alpha = .5, size = .8) +
scale_color_manual(values = c("FALSE" = NA,#248232
"TRUE" = "#D32934FF"),
labels = NULL) +
scale_fill_manual("", values = c("Favelas e lotes irregulares" = "#2BAA92FF")) +
labs(title = "<span style='font-size: 34pt;'>População em <span style = 'color:#D32934FF;'>áreas sem registro de IPTU</span> geralmente estão em <br><span style = 'color:#2BAA92FF;'>favelas ou lotes irregulares</span> (cada ponto representa 1000 pessoas)</span>") +
theme_void() +
theme(plot.title = element_markdown(), legend.position = "none")
# scale_colour_paletteer_d("lisa::AndyWarhol_2")
ggsave("tex/imagens/mapa_pontos.pdf", gg, width = 16, height = 20)
gg

Análise de resultados via join
df <- IPTU.censo %>%
filter(residencial == 1) %>%
st_drop_geometry() %>%
group_by(id_setor_censitario) %>%
summarize(# Setor censitário
populacao = median(v0001),
area_setor = median(area_setor) %>% as.numeric(),
domicilios = median(v0002),
domicilios_ocupados = median(v0007),
# Lote
unidades = sum(unidades * percent_intersec),
area_terreno = sum(area_terreno * percent_intersec),
area_construida = sum(area_construida * percent_intersec),
area_ocupada = sum(area_ocupada * percent_intersec),
) %>%
mutate(densidade = populacao * 10 ^ 6 / area_setor,
cota_parte = area_terreno / unidades,
CA = area_construida / area_terreno,
verticalizacao = area_construida / area_ocupada,
espectro_irregularidade = unidades / (unidades + domicilios))
modelsummary::datasummary_skim(df)
tinytable_sdx0eu8bx1fdzdk9eo65
| |
Unique |
Missing Pct. |
Mean |
SD |
Min |
Median |
Max |
Histogram |
| populacao |
1120 |
0 |
420.0 |
213.4 |
0.0 |
397.0 |
2221.0 |
 |
| area_setor |
18531 |
0 |
37267.5 |
109998.2 |
196.3 |
24549.0 |
6945494.4 |
 |
| domicilios |
543 |
0 |
189.6 |
88.1 |
0.0 |
182.0 |
1372.0 |
 |
| domicilios_ocupados |
471 |
0 |
162.5 |
75.3 |
0.0 |
156.0 |
950.0 |
 |
| unidades |
18524 |
0 |
141.1 |
97.1 |
0.0 |
130.1 |
1343.0 |
 |
| area_terreno |
18531 |
0 |
14702.9 |
18405.7 |
0.0 |
11595.9 |
1445351.6 |
 |
| area_construida |
18531 |
0 |
18094.3 |
13503.4 |
0.0 |
15739.1 |
163349.7 |
 |
| area_ocupada |
18531 |
0 |
7008.5 |
6150.2 |
0.0 |
5626.0 |
64738.0 |
 |
| densidade |
18341 |
0 |
28519.9 |
36697.1 |
0.0 |
17034.3 |
1162304.8 |
 |
| cota_parte |
18156 |
0 |
1282.3 |
14374.8 |
1.8 |
133.4 |
516201.6 |
 |
| CA |
18238 |
0 |
2.3 |
2.4 |
0.0 |
1.0 |
27.5 |
 |
| verticalizacao |
17941 |
0 |
4.7 |
5.5 |
1.0 |
1.8 |
77.4 |
 |
| espectro_irregularidade |
18360 |
0 |
0.4 |
0.2 |
0.0 |
0.4 |
1.0 |
 |
df %>%
select(id_setor_censitario, populacao, domicilios, unidades, area_setor) %>%
bind_rows(censo %>%
st_drop_geometry() %>%
select(id_setor_censitario, populacao = v0001, domicilios = v0002, area_setor) %>%
mutate(unidades = 0) %>%
anti_join(df %>% select(id_setor_censitario))) %>%
filter(populacao > 0) %>%
mutate(densidade = populacao/area_setor,
taxa_irregular = (domicilios - unidades) / domicilios,
espectro_irregularidade = unidades / (unidades + domicilios)) %>%
lm(log(densidade) ~ espectro_irregularidade, .) %>%
summary()
##
## Call:
## lm(formula = log(densidade) ~ espectro_irregularidade, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.0873 -0.4522 0.0812 0.7049 4.1956
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.10609 0.01168 -351.450 < 2e-16 ***
## espectro_irregularidade 0.14028 0.03279 4.278 1.89e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.224 on 26833 degrees of freedom
## Multiple R-squared: 0.0006816, Adjusted R-squared: 0.0006443
## F-statistic: 18.3 on 1 and 26833 DF, p-value: 1.893e-05
modelo1 <- df %>%
lm(densidade ~ cota_parte + verticalizacao + CA, data = .) %>%
summary()
modelo2 <- df %>%
filter(abs(espectro_irregularidade - .5) < .25) %>%
lm(densidade ~ cota_parte + verticalizacao + CA, data = .) %>%
summary()
modelo3 <- df %>%
filter(abs(espectro_irregularidade - .5) < .1) %>%
lm(densidade ~ cota_parte + verticalizacao + CA, data = .) %>%
summary()
modelsummary::modelsummary(
list("Todos setores" = modelo1, "Balanço entre 25% e 75%" = modelo2, "Balanço entre 40% e 60%" = modelo3),
gof_omit = "RMSE", estimate = "{estimate} {stars}",
statistic = "({std.error})")
tinytable_azsnqrglw69k1oc87i0x
| |
Todos setores |
Balanço entre 25% e 75% |
Balanço entre 40% e 60% |
| (Intercept) |
8603.264 *** |
4242.744 *** |
348.108 |
| |
(311.878) |
(347.853) |
(608.823) |
| cota_parte |
0.095 *** |
-0.018 |
-4.176 * |
| |
(0.016) |
(0.057) |
(1.712) |
| verticalizacao |
1058.806 *** |
1146.076 *** |
1283.914 *** |
| |
(59.059) |
(61.902) |
(76.170) |
| CA |
6524.131 *** |
7221.037 *** |
7874.627 *** |
| |
(132.933) |
(138.888) |
(171.647) |
| Num.Obs. |
18531 |
15616 |
9363 |
| R2 |
0.314 |
0.374 |
0.410 |
| R2 Adj. |
0.314 |
0.374 |
0.409 |
split <- rsample::initial_split(df, prop = 0.7)
train <- rsample::training(split)
test <- rsample::testing(split)
rf_model <- ranger::ranger(densidade ~ .,
data = train, num.trees = 5000, importance = "permutation")
## Growing trees.. Progress: 93%. Estimated remaining time: 2 seconds.
var.importance <- ranger::importance(rf_model)
gg <- ggplot(NULL, aes(x = var.importance,
y = fct_reorder(names(var.importance) %>% str_replace_all("_", " "),
var.importance))) +
geom_col(fill = "darkblue", alpha = 0.75) +
labs(x = "Importância", y = "", title = "Variáveis mais relevantes para explicar a densidade",
subtitle = "Resultados da metodologia via join") +
theme_classic()
ggsave("tex/imagens/importancia.pdf", gg, width = 7, height = 4)
remove(rf_model)
gg

Raster
geoms.IPTU <- IPTU.lote %>%
filter(! st_is_empty(geometry)) %>%
st_as_sf() %>%
mutate(area = st_area(geometry) %>% as.numeric()) %>%
select(id = setor, area, unidades, area_construida, area_ocupada, area_terreno)
geoms.censo <- censo %>%
mutate(area = st_area(geometry) %>% as.numeric()) %>%
select(id = id_setor_censitario, area, populacao = v0001, domicilios = v0002, domicilios_ocupados = v0007)
bbox <- st_bbox(geoms.censo)
raster <- raster::raster(xmn = bbox["xmin"], xmx = bbox["xmax"], ymn = bbox["ymin"], ymx = bbox["ymax"],
res = 800, crs = st_crs(geoms.censo)$proj4string, vals = NA)
result.IPTU <- exact_extract(raster, geoms.IPTU, include_xy = T,
include_cols = c("id", "area", "unidades", "area_construida", "area_ocupada",
"area_terreno"),
force_df = T, coverage_area = T, progress = F) %>%
bind_rows %>%
mutate(percent_intersect = coverage_area / area,
across(unidades:area_terreno, ~ .x * percent_intersect)) %>%
group_by(x, y) %>%
summarize(unidades = sum(unidades),
area_ocupada = sum(area_ocupada),
area_construida = sum(area_construida),
area_terreno = sum(area_terreno)) %>%
ungroup() %>%
mutate(cota_parte = area_terreno / unidades,
cota_parte_inversa = unidades / area_terreno,
CA = area_construida / area_terreno,
verticalizacao = area_construida / area_ocupada)
result.censo <- exact_extract(raster, geoms.censo, include_xy = T,
include_cols = c("id", "area", "populacao", "domicilios", "domicilios_ocupados"),
force_df = T, coverage_area = T, progress = F) %>%
bind_rows %>%
mutate(percent_intersect = coverage_area / area,
across(c(area, populacao, domicilios, domicilios_ocupados), ~ .x * percent_intersect)) %>%
group_by(x, y) %>%
summarize(populacao = sum(populacao),
area = sum(area),
domicilios = sum(domicilios),
domicilios_ocupados = sum(domicilios_ocupados))
result <- full_join(result.IPTU, result.censo) %>%
mutate(espectro_irregularidade = unidades / (unidades + domicilios),
densidade_residencial = populacao / area_terreno,
densidade_total = populacao / area)
result %>%
modelsummary::datasummary_skim()
tinytable_zbceqtd40kpw2vvpvl2e
| |
Unique |
Missing Pct. |
Mean |
SD |
Min |
Median |
Max |
Histogram |
| x |
59 |
0 |
331929.5 |
11037.7 |
313794.9 |
329794.9 |
360194.9 |
 |
| y |
91 |
0 |
7383596.5 |
18725.5 |
7343788.7 |
7388588.7 |
7415788.7 |
 |
| unidades |
1256 |
52 |
2083.5 |
2099.4 |
0.0 |
1733.9 |
21056.9 |
 |
| area_ocupada |
1256 |
52 |
103486.6 |
64038.7 |
1.5 |
113195.8 |
251213.6 |
 |
| area_construida |
1256 |
52 |
267192.3 |
267862.6 |
1.5 |
216295.9 |
2321166.4 |
 |
| area_terreno |
1256 |
52 |
217116.1 |
125745.5 |
7.8 |
238577.3 |
840996.4 |
 |
| cota_parte |
1254 |
52 |
2647.1 |
42935.6 |
6.5 |
140.7 |
1415766.7 |
 |
| cota_parte_inversa |
1255 |
52 |
0.0 |
0.0 |
0.0 |
0.0 |
0.2 |
 |
| CA |
1255 |
52 |
1.3 |
1.2 |
0.0 |
0.9 |
9.6 |
 |
| verticalizacao |
1246 |
52 |
2.6 |
2.3 |
1.0 |
1.8 |
28.9 |
 |
| populacao |
2403 |
0 |
4367.6 |
4804.2 |
0.0 |
2620.9 |
29598.3 |
 |
| area |
1787 |
0 |
580109.1 |
157947.2 |
0.1 |
640000.0 |
640000.0 |
 |
| domicilios |
2407 |
0 |
1905.6 |
2154.4 |
0.0 |
1108.2 |
17874.7 |
 |
| domicilios_ocupados |
2403 |
0 |
1646.2 |
1853.0 |
0.0 |
925.5 |
13732.5 |
 |
| espectro_irregularidade |
1256 |
52 |
0.4 |
0.2 |
0.0 |
0.4 |
1.0 |
 |
| densidade_residencial |
1256 |
52 |
0.6 |
9.2 |
0.0 |
0.0 |
271.1 |
 |
| densidade_total |
2354 |
0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
 |
gg <- result %>%
pivot_longer(c(CA, cota_parte, verticalizacao, densidade_total)) %>%
group_by(name) %>%
mutate(value = ntile(value, 10) %>% as.factor(),
name = case_when(name == "CA" ~ "Densidade Construtiva",
name == "cota_parte" ~ "Densidade Habitacional",
name == "densidade_total" ~ "Densidade populacional",
name == "verticalizacao" ~ "Verticalização",
.default = NA)) %>%
ggplot() +
geom_sf(data = distrito) +
geom_tile(aes(x = x, y = y, fill = value), alpha = .75, color = NA, lwd = .01) +
# facet_wrap(~name, ncol = 4) +
theme_void() +
theme(strip.text = element_text(size = 20)) +
labs(fill = "Decil") +
scale_fill_viridis_d()
ggsave("tex/imagens/rasters.pdf", gg + facet_wrap(~name, ncol = 2), width = 8, height = 12)
ggsave("tex/imagens/rasters_wide.pdf", gg + facet_wrap(~name, ncol = 4), width = 16, height = 6)
gg + facet_wrap(~name, ncol = 2)

gg <- result |>
mutate(corte = "todos") |>
bind_rows(result |> filter(abs(espectro_irregularidade - .5) <= .4) |> mutate(corte = "<.4")) |>
bind_rows(result |> filter(abs(espectro_irregularidade - .5) <= .3) |> mutate(corte = "<.3")) |>
bind_rows(result |> filter(abs(espectro_irregularidade - .5) <= .2) |> mutate(corte = "<.2")) |>
bind_rows(result |> filter(abs(espectro_irregularidade - .5) <= .1) |> mutate(corte = "<.1")) |>
bind_rows(result |> filter(abs(espectro_irregularidade - .5) <= .05) |> mutate(corte = "<.05")) |>
filter(area > max(area) * .99) |>
mutate(value = ntile(densidade_residencial, 10) %>% as.factor()) |>
ggplot() +
geom_sf(data = distrito) +
geom_tile(aes(x = x, y = y, fill = value), alpha = .75, color = NA, lwd = .01) +
facet_wrap(~corte) +
theme_void() +
theme(strip.text = element_text(size = 12)) +
labs(fill = "Decil") +
scale_fill_viridis_d()
ggsave("tex/imagens/rasters_densidade.pdf", gg, width = 10, height = 10)
gg

gg <- result %>%
ggplot(aes(x = x, y = y)) +
geom_raster(aes(fill = populacao)) +
theme(
legend.position = "none",
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank()
) +
labs(title = "População em São Paulo (2022)", x = "", y = "") +
coord_equal() +
scale_fill_viridis_c()
# plot_gg(gg, width = 5, height = 6, raytrace = FALSE, preview = TRUE)
plot_gg(gg, multicore = TRUE, scale = 150, width = 5, height = 7, pointcontract = 1,
theta = -50, phi = 40, zoom = .65, fov = 8, windowsize = c(3840, 2160))
render_snapshot("tex/imagens/mapa3d", software_render = TRUE, width = 3840, height = 2160,
camera_lookat = c(50, -125, 50))
(result %>%
filter(area > max(area) * .99) %>%
arrange(desc(densidade_total)) %>%
mutate(top = ifelse(row_number() <= 5, row_number(), NA) %>% as.factor()) %>%
ggplot() +
geom_sf(data = distrito, color = "grey", aes(text = ds_nome)) +
geom_tile(aes(x = x, y = y, fill = top, text = paste("pop: ", populacao, "\n",
"unidades: ", unidades, "\n",
"domicilios: ", domicilios, "\n",
"densidade: ", densidade_total, "\n",
"area: ", area,
sep = ""))) +
scale_fill_viridis_d() +
theme_void()) %>%
plotly::ggplotly()
tabela <- result %>%
filter(area > max(area) * .99) %>%
mutate(across(everything(), ~ replace_na(., 0)),
across(c(densidade_residencial, densidade_total), ~ . * 10^6)) %>%
arrange(desc(densidade_total)) %>%
head(5) %>%
mutate(localizacao = c("1. Paraisópolis", "2. Heliópolis ", "3. Sé (Bela Vista)", "4. Paraisópolis", "5. Heliópolis")) %>%
select("Localização" = localizacao, "População" = populacao, "Domicílios (Censo)" = domicilios,
"Domicílios Ocupados" = domicilios_ocupados, "Unidades (IPTU)" = unidades,
"Espectro irregularidade" = espectro_irregularidade, "Densidade populacional" = densidade_total, "Área" = area) %>%
pivot_longer(cols = -Localização, names_to = "Variável") %>%
pivot_wider(id_cols = Variável, names_from = Localização, values_from = value) %>%
gt() %>%
tab_spanner("Favelas", columns = c(2,3,5,6)) %>%
fmt_number(rows = c(1,2,3, 4, 6, 7), decimals = 0, sep_mark = ".", dec_mark = ",") %>%
fmt_percent(rows = 5)
# gtsave(tabela %>% tab_options(table.font.size = 12), "tex/tabelas/top5.tex")
# workaround bug da gt quando exporta para latex
as_latex(tabela %>% tab_options(table.font.size = 13)) |>
as.character() |>
sub(pattern = "(\\\\begin\\{longtable\\}\\{lrrrrr\\})",
replacement = "\\1 \n\\\\caption{Células do \\\\textit{raster} que apresentam a maior densidade populacional}\n\\\\label{tab:top5}\\\\\\\\",
x = _) |>
write_lines("tex/tabelas/top5.tex")
gtsave(tabela %>% tab_options(table.font.size = 9), "tex/tabelas/top5_small.tex")
tabela
| Variável |
Favelas
|
3. Sé (Bela Vista) |
| 1. Paraisópolis |
2. Heliópolis |
4. Paraisópolis |
5. Heliópolis |
| População |
29.598 |
25.280 |
23.824 |
22.920 |
24.576 |
| Domicílios (Censo) |
11.655 |
10.178 |
9.361 |
9.001 |
17.875 |
| Domicílios Ocupados |
10.791 |
9.269 |
8.810 |
8.274 |
13.732 |
| Unidades (IPTU) |
0 |
1.857 |
7 |
3 |
21.057 |
| Espectro irregularidade |
0.00% |
15.43% |
0.08% |
0.03% |
54.09% |
| Densidade populacional |
46.247 |
39.500 |
37.225 |
35.813 |
38.400 |
| Área |
640.000 |
640.000 |
640.000 |
640.000 |
640.000 |
gg <- result %>%
drop_na() %>%
ggplot() +
geom_sf(data = distrito, color = NA) +
geom_tile(aes(x = x, y = y, fill = espectro_irregularidade), alpha = .75, color = "grey") +
theme_void() +
scale_fill_gradient2(high = "#2F191B", low = "#D32934", mid = "white", na.value = NA,
midpoint = .5, limits = c(0,1)) +
theme(legend.position = c(.8, .35), legend.title = element_markdown(size = 25),
plot.title = element_markdown(), legend.text = element_text(size = 25),
legend.key.size = unit(1.75, 'cm')) +
labs(fill = "Espectro de irregularidade", title = "<span style='font-size: 35pt;'>Áreas com <span style = 'color:#D32934;'>menos domicílios registrados no IPTU, comparados <br>ao censo</span> geralmente estão em regiões afastadas do centro</span>")
ggsave("tex/imagens/balanco_raster.pdf", gg, width = 16, height = 20)
gg

gg <- ggcorrplot::ggcorrplot(result |>
filter(abs(espectro_irregularidade - .5) < .1) |>
mutate(tx_ocupacao = area_ocupada / area_terreno) |>
select("Taxa de ocupação" = tx_ocupacao,
"Verticalização" = verticalizacao,
"Densidade Habitacional" = cota_parte,
"Densidade Construtiva" = CA,
"Densidade Total" = densidade_total,
"Densidade Residencial" = densidade_residencial) |>
drop_na() |>
cor(),
show.diag = FALSE, type = "upper", legend.title = "Correlação",
colors = c("red", "white", "blue"), lab = TRUE, outline.color = "white") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
ggsave("tex/imagens/corrplot.pdf", gg, width = 7, height = 4)
gg

Análise resultados via Raster
modelo1 <- result %>%
mutate(across(c(densidade_residencial, densidade_total), ~ . * 10^6)) %>%
filter(abs(espectro_irregularidade - .5) < .15, populacao > 0) %>%
lm(densidade_residencial ~ cota_parte + verticalizacao + CA, data = .) %>%
summary()
modelo2 <- result %>%
mutate(across(c(densidade_residencial, densidade_total), ~ . * 10^6)) %>%
filter(abs(espectro_irregularidade - .5) < .15, populacao > 0) %>%
lm(log(densidade_residencial) ~ cota_parte + verticalizacao + CA, data = .) %>%
summary()
modelo3 <- result %>%
mutate(across(c(densidade_residencial, densidade_total), ~ . * 10^6)) %>%
filter(abs(espectro_irregularidade - .5) < .05) %>%
lm(densidade_residencial ~ cota_parte + verticalizacao + CA, data = .) %>%
summary()
modelo4 <- result %>%
mutate(across(c(densidade_residencial, densidade_total), ~ . * 10^6)) %>%
filter(abs(espectro_irregularidade - .5) < .05,
populacao > 0) %>%
lm(log(densidade_residencial) ~ cota_parte + verticalizacao + CA, data = .) %>%
summary()
modelos <- list("Nível (A)" = modelo1,
"Log (B)" = modelo2,
"Nível (C)" = modelo3,
"Log (D)" = modelo4)
tabela <- modelsummary(modelos,
estimate = c("Coeficiente" = "estimate"),
statistic = c("*" = "stars"),
shape = term ~ model + statistic,
output = "gt",
align = "lrlrlrlrl",
coef_map = c("cota_parte" = "Densidade Habitacional",
"CA" = "Densidade Construtiva",
"verticalizacao" = "Verticalização"),
add_rows = modelos %>%
sapply(get_gof) %>%
rbind %>%
as_tibble(rownames = 'variable') %>%
unnest(cols = c(`Nível (A)`, `Log (B)`, `Nível (C)`, `Log (D)`)) %>%
mutate(across(c(`Nível (A)`, `Log (B)`, `Nível (C)`, `Log (D)`), ~
ifelse(variable == "nobs", round(.) %>% as.character(), sprintf('%.3f',.))),
variable = case_when(variable == "r.squared" ~ "R2",
variable == "adj.r.squared" ~ "R2 ajustado",
variable == "nobs" ~ "Observações",
.default = NA)) %>%
drop_na() %>%
mutate(!!!list("As" = "", "Bs" = "", "Cs" = "", "Ds" = "")) %>%
select(variable, "Nível (A)", As, "Log (B)", Bs, "Nível (C)", Cs, "Log (D)", Ds,)) %>%
tab_spanner("Espectro irregularidade: 35 a 65%", columns = 2:5) %>%
tab_spanner("Espectro irregularidade: 45 a 55%", columns = 6:9) %>%
tab_style(style = cell_borders(sides = "top", weight = px(0.5)),
locations = cells_body(rows = 5)) %>%
sub_values(values = "+", replacement = ".")
gtsave(tabela %>% tab_options(table.font.size = 12), "tex/tabelas/regressao.tex")
gtsave(tabela %>% tab_options(table.font.size = 10), "tex/tabelas/regressao_small.tex")
tabela
|
Espectro irregularidade: 35 a 65%
|
Espectro irregularidade: 45 a 55%
|
| |
Nível (A)
|
Log (B)
|
Nível (C)
|
Log (D)
|
| Coeficiente |
* |
Coeficiente |
* |
Coeficiente |
* |
Coeficiente |
* |
| Densidade Habitacional |
-9.892 |
*** |
-0.002 |
*** |
-9.410 |
** |
-0.002 |
*** |
| Densidade Construtiva |
10951.872 |
*** |
0.197 |
*** |
11470.672 |
*** |
0.274 |
*** |
| Verticalização |
816.466 |
. |
0.014 |
|
1398.701 |
** |
0.024 |
. |
| R2 |
0.507 |
|
0.638 |
|
0.639 |
|
0.737 |
|
| R2 ajustado |
0.505 |
|
0.636 |
|
0.636 |
|
0.735 |
|
| Observações |
741 |
|
741 |
|
332 |
|
332 |
|
pvals <- tibble()
corte <- .05
for (corte in seq(.005, .5, length.out = 1000)){
pval.nivel <- result |>
filter(abs(espectro_irregularidade - .5) < corte) |>
lm(densidade_residencial ~ cota_parte + verticalizacao + CA, data = _) |>
summary()
pval.log <- result |>
filter(abs(espectro_irregularidade - .5) < corte) |>
lm(log(densidade_residencial) ~ cota_parte + verticalizacao + CA, data = _) |>
summary()
pvals <- tibble(corte = corte,
nobs = pval.nivel$df[2] + 4,
cota_parte = pval.nivel$coefficients[,4][["cota_parte"]],
verticalizacao = pval.nivel$coefficients[,4][["verticalizacao"]],
CA = pval.nivel$coefficients[,4][["CA"]],
forma = "nivel") |>
bind_rows(pvals)
pvals <- tibble(corte = corte,
nobs = pval.log$df[2] + 4,
cota_parte = pval.log$coefficients[,4][["cota_parte"]],
verticalizacao = pval.log$coefficients[,4][["verticalizacao"]],
CA = pval.log$coefficients[,4][["CA"]],
forma = "log") |>
bind_rows(pvals)
}
gg <- pvals |>
select(corte, forma,
"Densidade Contrutiva" = CA,
"Densidade Habitacional" = cota_parte,
"Verticalização" = verticalizacao) |>
pivot_longer(3:5) |>
ggplot(aes(x = corte, y = value, color = name)) +
geom_hline(yintercept = 0) +
geom_line(lwd = 1, alpha = .8) +
geom_hline(yintercept = .01, linetype = "dashed") +
facet_grid(rows = "forma") +
theme_classic() +
scale_y_sqrt(labels = scales::percent, breaks = c(.01, .05, .1, .5, 1), expand = c(0, 0)) +
scale_x_continuous(labels = scales::percent, expand = c(0, 0)) +
labs(x = "Corte do espectro de irregularidade",
y = "p.valor (escala em raiz quadrada)",
colour = "Variável")
ggsave("tex/imagens/pvals.pdf", gg, width = 8, height = 6)
gg

gg <- pvals |>
filter(forma == "log") |>
ggplot(aes(x = corte, y = nobs)) +
geom_line() +
geom_hline(yintercept = min(pvals$nobs), linetype = "dashed", alpha = .5) +
theme_classic() +
scale_y_continuous(limits = c(0,1250), breaks = c(pretty(pvals$nobs), min(pvals$nobs))) +
scale_x_continuous(labels = scales::percent) +
labs(x = "Corte do espectro de irregularidade", y = "Número de observações")
ggsave("tex/imagens/nobs.pdf", gg, width = 6, height = 4)
gg

split <- rsample::initial_split(result %>%
drop_na() %>%
filter(abs(espectro_irregularidade - .5) <.1,
area > max(area) * .99) %>%
select(-c(densidade_total, populacao, area,
cota_parte_inversa, x, y, espectro_irregularidade)),
prop = 0.7)
train <- rsample::training(split)
test <- rsample::testing(split)
rf_model <- ranger::ranger(densidade_residencial ~ .,
data = train, num.trees = 100000, importance = "permutation")
rf_preds <- predict(rf_model, data = test)$predictions
rf_model_restrito <- ranger::ranger(densidade_residencial ~ CA + cota_parte + verticalizacao,
data = train, num.trees = 100000, importance = "permutation")
rf_preds_restrito <- predict(rf_model_restrito, data = test)$predictions
gg <- tibble(names = names(rf_model_restrito$variable.importance),
importance = rf_model_restrito$variable.importance) |>
mutate(names = case_when(names == "CA" ~ "Densidade Construtiva",
names == "cota_parte" ~ "Densidade Habitacional",
names == "verticalizacao" ~ "Verticalização",
.default = names)) |>
ggplot(aes(x = importance,
y = fct_reorder(names, importance))) +
geom_col(fill = "darkblue", alpha = 0.75) +
labs(x = "Importância", y = "Variável", title = "") +
theme_classic()
ggsave("tex/imagens/var_importance_restrito.pdf", gg, width = 5.5, heigh = 2)
gg

gg <- tibble(names = names(rf_model$variable.importance),
importance = rf_model$variable.importance) |>
mutate(names = case_when(names == "CA" ~ "Densidade Construtiva",
names == "cota_parte" ~ "Densidade Habitacional",
names == "verticalizacao" ~ "Verticalização",
names == "area_terreno" ~ "Área Terreno",
names == "area_ocupada" ~ "Área Ocupada",
names == "area_construida" ~ "Área Construída",
names == "unidades" ~ "Unidades",
names == "domicilios" ~ "Domicílios",
names == "domicilios_ocupados" ~ "Domicílios Ocupados",
.default = names)) |>
ggplot(aes(x = importance,
y = fct_reorder(names, importance))) +
geom_col(fill = "darkblue", alpha = 0.75) +
labs(x = "Importância", y = "Variável", title = "") +
theme_classic()
ggsave("tex/imagens/var_importance.pdf", gg, width = 7, heigh = 4)
gg

df <- tibble()
for (cota_parte in seq(5, 100, by = 5)){
for (CA in seq(0, 8, by = .5)){
for (verticalizacao in c(8, 20)){
df <- df |>
bind_rows(tibble(cota_parte = cota_parte, CA = CA, verticalizacao = verticalizacao))
}
}
}
gg <- df |>
mutate(preds = predict(rf_model_restrito,
data = tibble(CA = CA,
cota_parte = cota_parte,
verticalizacao = verticalizacao)) |>
ranger::predictions(),
preds = cut(preds * 5000, breaks = c(0:10) * 100, dig.lab = 6),
verticalizacao = paste(verticalizacao, "pavimentos") |>
factor(levels = c("8 pavimentos", "20 pavimentos"))) |>
ggplot(aes(x = CA, y = cota_parte, fill = preds)) +
geom_tile() +
facet_wrap(~verticalizacao, dir = "h") +
scale_fill_viridis_d() +
theme_classic() +
labs(fill = "População Prevista", x = "Coeficiente de Aproveitamento (CA)",
y = "Cota Parte")
ggsave("tex/imagens/previsoes.pdf", width = 9, height = 4.25)
gg

split <- rsample::initial_split(result |>
filter(abs(espectro_irregularidade - .5) < 1) |>
select(densidade_total,
dens_cons = CA,
dens_hab = cota_parte,
vertic = verticalizacao), prop = 0.7)
rtree <- rpart::rpart(densidade_total ~ dens_cons + dens_hab + vertic,
data = rsample::training(split))
visNetwork::visTree(rtree, direction = "LR", collapse = TRUE, legend = FALSE)
lm_model <- lm(densidade_residencial ~ CA + cota_parte + verticalizacao, data = train)
lm_preds <- predict(lm_model, newdata = test)
1 - sum((rf_preds - test$densidade_residencial)^2)/sum((test$densidade_residencial - mean(test$densidade_residencial))^2)
## [1] 0.8473629
1 - sum((rf_preds_restrito - test$densidade_residencial)^2)/sum((test$densidade_residencial - mean(test$densidade_residencial))^2)
## [1] 0.8196794
1 - sum((lm_preds - test$densidade_residencial)^2)/sum((test$densidade_residencial - mean(test$densidade_residencial))^2)
## [1] 0.5063794
split <- rsample::initial_split(result %>%
drop_na() %>%
filter(abs(espectro_irregularidade - .5) <.1,
area > max(area) * .99) %>%
select(-c(densidade_residencial, populacao, area,
cota_parte_inversa, x, y, espectro_irregularidade)),
prop = 0.7)
train <- rsample::training(split)
test <- rsample::testing(split)
rf_model <- ranger::ranger(densidade_total ~ .,
data = train, num.trees = 100000, importance = "permutation")
rf_preds <- predict(rf_model, data = test)$predictions
rf_model_restrito <- ranger::ranger(densidade_total ~ CA + cota_parte + verticalizacao,
data = train, num.trees = 100000, importance = "permutation")
rf_preds_restrito <- predict(rf_model_restrito, data = test)$predictions
gg <- tibble(names = names(rf_model$variable.importance),
importance = rf_model$variable.importance) |>
mutate(names = case_when(names == "CA" ~ "Densidade Construtiva",
names == "cota_parte" ~ "Densidade Habitacional",
names == "verticalizacao" ~ "Verticalização",
names == "area_terreno" ~ "Área Terreno",
names == "area_ocupada" ~ "Área Ocupada",
names == "area_construida" ~ "Área Construída",
names == "unidades" ~ "Unidades",
names == "domicilios" ~ "Domicílios",
names == "domicilios_ocupados" ~ "Domicílios Ocupados",
.default = names)) |>
ggplot(aes(x = importance,
y = fct_reorder(names, importance))) +
geom_col(fill = "darkblue", alpha = 0.75) +
labs(x = "Importância", y = "Variável", title = "") +
theme_classic()
ggsave("tex/imagens/var_importance_densmod.pdf", gg, width = 7, heigh = 4)
gg

lm_model <- lm(densidade_total ~ CA + cota_parte + verticalizacao, data = train)
lm_preds <- predict(lm_model, newdata = test)
1 - sum((rf_preds - test$densidade_total)^2)/sum((test$densidade_total - mean(test$densidade_total))^2)
## [1] 0.9618918
1 - sum((rf_preds_restrito - test$densidade_total)^2)/sum((test$densidade_total - mean(test$densidade_total))^2)
## [1] 0.4427002
1 - sum((lm_preds - test$densidade_total)^2)/sum((test$densidade_total - mean(test$densidade_total))^2)
## [1] 0.2868002